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Abstract 


Trailing vortices generated by lifting surfaces such as helicopter rotor blades, 
ship propellers, fixed wings, and canard control surfaces are known to be the 
source of noise, vibration, cavitation, degradation of performance, and other 
hazardous problems. Controlling these vortices is therefore of practical interest. 
The formation and behavior of the trailing vortices are studied in the present 
research. In addition, wing-tip blowing concepts employing axial blowing and 
spanwise blowing are studied to determine their effectiveness in controlling these 
vortices and their effects on the performance of the wing. The three- 
dimensional, unsteady, thin-layer compressible Navier-Stokes equations are 
solved using a time-accurate, implicit, finite difference scheme that employs LU- 
ADI factorization. The wing-tip blowing is simulated using the actuator plane 
concept, thereby not requiring resolution of the jet slot geometry. Furthermore, 
the solution blanking feature of the chimera scheme is used to simplify the 
parametric study procedure for the wing-tip blowing. Computed results are 
shown to compare favorably with experimental measurements. It is found that 
axial wing-tip blowing, although delaying the rolling-up of the trailing vortices 
and the near-field behavior of the flowfield, does not dissipate the circulation 
strength of the trailing vortex farther downstream. Spanwise wing-tip blowing 




has the effect of displacing the trailing vortices outboard and upward. The 
increased “wing-span” due to the spanwise wing-tip blowing has the effect of lift 
augmentation on the wing and the strengthening of the trailing vortices. 
Secondary trailing vortices are created at high spanwise wing-tip blowing 


intensities. 
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Chapter 1 


Introduction 


1.1 Motivation 

Near the tip of a lifting surface, the pressure differential between the upper and 
lower surfaces causes the fluid to move around the tip from the higher pressure 
lower surface towards the lower pressure upper surface. As the fluid undergoes 
this highly 3-dimensional maneuver at the tip, it realigns the bound vortex 
(caused by the boundary layer vorticity of the lifting surface) which is parallel to 
the span of the lifting surface into the downstream direction. This tip vortex is 
convected downstream and is typically called the trailing vortex. 

Trailing vortices of a conventional fixed-wing aircraft are depicted in 
Figure 1.1. The vortices generated by larger fixed-wing aircraft may cause 
undesirable and sometimes uncontrollable rolling moment and other hazardous 
effects on smaller aircraft flying into their wakes. (See Figure 1.2.) These 
trailing vortices may persist for a considerable distance; and result in a rather 
long safety distance (thus long interval) required between flights. The problem 
can be severe especially during take-off and landing at a busy airport. Trailing 
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vortices from the canards of many modem fighter aircraft can also pose 
problems on the wings of the aircraft resulting in the degradation of 
performance. 

For helicopter rotor blades, trailing vortices are convected downward in a 
helical manner (see Figure 1 .3) due to the downwash and the rotational motion 
of the lifting surfaces. Vortices shed from the preceding blades of a helicopter 
may interact with the following blades (called blade-vortex interaction; see 
Figure 1.4) during forward and descending flights. This is a major source of 
noise and vibration for helicopters (see, for example, Schmitz and Yu [1983]). 
Ship propellers have a very similar problem encountered by helicopter rotor 
blades. In addition, erosion on propeller blades may occur as the propeller 
blades strike the low pressure air pockets (called cavitation) created by the 
trailing vortices (see Figure 1.5). 

It is therefore of practical interest to explore the effects of control devices 
such as axial and spanwise wing-tip blowing on the formation and behavior of tip 
vortices and on the performance of the wing. Spanwise blowing has been studied 
also as a means of control to effect a change of vortex position on one wing-tip 
and thereby produce a rolling moment through asymmetric wing lift (see Tavella 
etal [1988]). 

For several reasons, then, it is desirable to better understand 
computationally the origin of the tip vortex and the wing it is influenced by 
blowing - i.e. whether it can be displaced to produce a rolling moment and 
whether its behavior downstream can be affected to reduce noise and reduce the 
hazard to following aircraft. 
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1.2 Previous Work 

1.2.1 Experimental and Theoretical Work 

The phenomenon of the trailing vortex from a lifting surface was recognized as 
early as 1907 by Lanchester [1907] (see Figure 1.6). Some subsequent landmark 
works include Westwater’s calculation of the vortex roll-up and trajectory using 
the two-dimensional (Trefftz plane) model of a vortex sheet consisting of point, 
doubly infinite, vortex elements (see Westwater [1935]). Betz [1932] analyzed 
the movement of the center of gravity of the trailing vortices in the two- 
dimensional cross-plane. Spreiter and Sacks [1951] studied the inviscid effects of 
lifting surface parameters such as aspect ratio, lift coefficient and wing chord on 
the downwash by solving the Laplace equation. 

By correlating with experimental data, McCormick et al [1968], Roberts 
[1975, 1984], Phillips [1981] proposed analytical scaling laws to describe the 
structure of viscous turbulent trailing vortices with some success. Various 
experimental studies involving scaled model testing using wind-tunnels and full- 
scale flight testing were also carried out especially in the 1970s (see, for 
example, Olsen et al [1971], and Hallock and Eberle [1977]). A more recent and 
comprehensive experimental study was performed by McAlister and Takahashi 
[1991]. 

Concurrently, efforts were being made to minimize the effects of trailing 
vortices. Most of the efforts employed various wing-tip configurations to 
spread, relocate, and dissipate the trailing vortices (see Figure 1.7). For 
example, the configuration with gulled outer panels in Figure 1.7 serves to 
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divide the total concentrated vorticity into two vortices, each of lesser strength, 
at two different spanwise locations, (one at the juncture of the gulled tip and 
main lifting surface and the other at the tip), and also to relocate the tip vortex 
with respect to a dimension normal to the wing chord plane because of the 
deflection of the tip of the gulled section (see White [1973]). The configuration 
with a drag spoiler fixed at the wing-tip is known to increase the aircraft drag 
(see Corsiglia et al [1971]). These fixed configurations, although they may be 
optimized for a particular flight condition, suffer degradation in performance 
for other flight conditions. On the other hand, configurations employing active 
control devices such as wing-tip blowing allow adjustments to be made as the 
flight condition changes by changing the blowing intensity. The emphasis of the 
present research effort is on these active control configurations employing wing- 
tip blowing. 

Some experimental investigations had been carried out to study the effects 
of axial wing-tip blowing (or mass injection) on the behavior of trailing vortices 
at low Reynolds number (see, for example, Snedeker [1971], Mason and 
Marchman [1973], White [1973], and Dunham [1976]). However, it was not 
conclusive whether this concept of employing axial wing-tip blowing was 
effective since conflicting observations were reported. 

The concept of employing spanwise wing-tip blowing had also been 
carried out experimentally (see, for example, Ayers and Wilde [1956], Wu et al 
[1984], Tavella et al [1985], and Lee et al [1989]). It was observed that lift 
augmentation could be achieved via spanwise wing-tip blowing. 
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Olsen et al [1971], Hallock and Eberle [1977], and Hallock [1992] provide 
more details on various research efforts on the trailing vortex problem. 

1.2.2 Computational Work 

A preliminary study on the tip vortex from a low aspect ratio wing was studied 
by Mansour [1985] using the Computational Fluid Dynamics (CFD) 
methodology. Srinivasan et al [1988] extended the study to include the effects of 
wing planform shape and the wing-tip shape on the structure of the trailing 
vortices. Strawn [1991] applied the unstructured adaptive-grid scheme to study 
the structure of inviscid trailing vortices by solving the Euler equations. Dacles- 
Mariani et al [1993] incorporated the experimental data of Chow et al [1993] into 
the computations to study the near-field structure of the turbulent wing-tip 
vortex. Childs [1986] did a preliminary study on the lift augmentation due to the 
spanwise wing-tip blowing. Wong and Kandil [1992] studied the effects of 
angled wing-tip blowing on the behavior of trailing vortices. 

1.3 Current Approach 

Despite the fact that a tremendous amount of work has been done on the trailing 
vortex problem, this complex flow phenomenon is still not very well understood, 
especially when active control devices such as wing-tip blowing are employed. 

The present study attempts to apply the CFD methodology to better 
understand the formation and behavior of trailing vortices and to evaluate the 
effects of axial and spanwise wing-tip blowing on these vortices for both fixed- 
wing and rotary-wing aircraft applications. However, no numerical 
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computations are performed for the rotorcraft configuration because of the 
enormous amount of computational effort required (see Srinivasan et al [1992]). 
A geometrically simple and computationally cheaper fixed-wing configuration is 
used instead. This is motivated by the findings of Tung et al [1983], who had 
investigated the structure of trailing vortices generated by model rotor blades 
and observed that the structure of the trailing vortices generated by rotor blades 

;Nr 

was very similar to that generated by fixed-wing aircraft, and Spivey [1968], 
who observed that the centrifugal effect due to the rotating rotor blade had very 
little impact on the path of the trailing vortices. 

1.4 Thesis Outline 

This thesis is divided into six chapters. Following the introduction in this 
chapter. Chapter 2 describes the equations governing the present research 
problem, the approximations introduced, the turbulence models used, the 
numerical algorithm and the artificial dissipation applied to solve these equations. 
Chapter 3 covers the grid generation method and grid topology used. Chapter 4 
discusses the implementation of the boundary conditions required by the Navier- 
Stokes equations for cases with and without wing-tip blowing. Chapter 5 
presents the computed results for the unblown cases, the axial and spanwise 
blowing cases and compares them with the available experimental data and 
theoretical analyses. The final chapter, Chapter 6, concludes the findings of this 
research and outlines possible improvements and future work. 
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Figure 1.1: Trailing vortices of a conventional fixed-wing aircraft (taken from 
Kuethe and Chow [1986]). 
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STRUCTURAL LOAD FACTORS 


Figure 1 .2: Hazardous effects on smaller aircraft flying into the wakes of larger 
aircraft (taken from Gessow [1976]). 
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Figure 1.3: Schematic sketch of the wake structure for a single rotor blade in 
hover (taken from Stepniewski and Keys [1984]). 
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Figure 1.4: Blade-vortex interaction of a helicopter during forward and/or 
descending flights (taken from White [1973]). 
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Figure 1.5: Progression of cavitation around propeller blades (taken from 

Trevena [1987]). 
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Figure 1.6: Lanchester’s drawing of the vortex roll-up (taken from Hackett and 
Evans [1971]). 
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Figure 1.7: Wing-tip configurations for controlling trailing vortices (taken from 
White [1973]). 




Chapter 2 


Governing Equations and 
Numerical Method 

2.1 Governing Equations and Approximations 

2.1.1 The Navier-Stokes Equations 

The governing equations used in the present trailing vortex study are the three- 
dimensional Navier-Stokes equations. Under the assumptions of no body forces 
and no external heat addition, the Navier-Stokes equations can be written in 
conservation-law form in nondimensional Cartesian coordinates as follows (see, 
for example, Peyret and Viviand [1975]): 

+ + + = + + ( 2 . 1 ) 

dt dx dy dz dx dy dz 

where the conservative flow variables vector Q is defined as 
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Q = 


pu 

pv 

pw 

e 


The inviscid or “Euler” flux vectors E, F, G are given by 


pu 


pv 


pw 

pu 2 +p 


puv 


puw 

puv 

, F = 

pv 2 +p 

i G = 

pvw 

puw 


pvw 


pw 2 + p 

_u(e + p)_ 


_v(e + p)_ 


w(e + p ) 


and the viscous flux vectors E v , F v , G„ are given by 


' 0 ' 


" 0 “ 


‘O' 
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ii 
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A. 


.ft. 


A. 


The components of the viscous stresses are 


* ia =Mu x + v y + w z ) + 2pu x 

ft = — ft*. + + VT^ + wr xt 

*yy=A'A + V y +W z )+2pV y 

T a = l(u x + v y +w z ) + 2pw t 

P> = jt 9 > e ‘ +ur y* + vr yy + wr y* 

r xy = r yz =p(u y +v z ) 

^ =//(^ + wJ 
T w = T w =/i(v z +W y ) 

ft - ^ ft*, + MT, + VT^ + 


( 2 . 2 ) 


( 2 . 3 ) 


( 2 . 4 ) 


( 2 . 5 ) 
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where u x = du/dx, v y = dv/dy, d x e t = dejdx, etc. 

The density p is nondimensionalized by the freestream density p„, the 
velocity components u, v and w by the freestream speed of sound a„, and the 
total energy per unit volume e by p„al. Conforming to the usual convention, u 
is in the wing chord (x) direction (positive aft), v is in the wing spanwise (y) 
direction (positive outboard), and w is in the vertical (z) direction (positive 
upwards). The coefficient of viscosity p is normalized by its freestream value, 
p_, and the time is normalized by c/a„ where c is the wing chord. The 
Reynolds number is defined as Re = p„,ac/p„ . The Prandtl number Pr is 
defined as Pr = c p p/k where c p is the specific heat at constant pressure and k is 

the coefficient of thermal conductivity. Also, y is the ratio of specific heats 
which for air is equal to 1 .4. 

The internal energy per unit mass e l is related to the total energy e as 
follows: 

e + ( 2.6) 


Pressure is related to the conservative flow variables through the equation of 
state for a perfect gas 


p = (y-l) 


-!<-* 


+ v 2 +w 2 ) 


(2.7) 


The fluid in the present study is also considered to be Newtonian (i.e., viscous 
stresses are linearly related to the rates of strain), isotropic (i.e., having no 
preferred direction), and that Stokes’ hypothesis, which states that the bulk 
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viscosity (A +%(i) is zero, is satisfied. These assumptions are found to be valid 
for aerodynamic applications. 


2.1.2 Coordinate Transformation 


Next, the Navier-Stokes equations are transformed from the Cartesian 
coordinates (x, y, z, t ) to a generalized, body-fitted, curvilinear coordinate 
system (£, rj, C,, t). This makes the formulation independent of the body 
geometry thereby easing the specification of the boundary conditions. The 
transformation from the physical domain to the computational domain also 
allows standard differencing schemes for equi-spaced grid points to be used for 
spatial derivatives. In addition, it also allows the thin-layer approximation (to be 
discussed in Section 2.1.3) to be applied in a straight-forward manner. The 
coordinate transformation is defined by: 


I = %(x,y,z,t) 
7] = rj(x,y,z,t) 
£ = £(x,y,z,t) 
z = t 


( 2 . 8 ) 


where t and r are independent variables of time in the physical and transformed 
coordinates, respectively. The airfoil surface in the chordwise direction is 
transformed to the £ -coordinate, the spanwise direction is transformed to the 77 - 
coordinate and the ^-coordinate is normal to the wing surface. Chain rule is 
applied in the transformation procedure which can be found in Viviand [1974] 
and Vinokur [1974]. The resulting transformation Jacobian is 
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/ = 1 / det 


l 

y* 

z. 


X i x ; 

y, y; 


(2.9) 


where x = dx/d % , x n = dx/drj, etc, and the transformation metrics are given by 


£ x =J(y n z ; -y ; z„) 

Z y =J( X rZ„- X r ,Z;) 


n* = * -y&) 

T iy = J ( x fr- x &) 


C, = ^(y { z n~y^ 

Zy = J(x tl Zs-x i z n ) 


~ X T ^ x y Z^y Z t £ z 
TJ l =-X T TJ x -y r TJy~Z T TJ z 

Z t = -x£*-y£ y -z£z 


( 2 . 10 ) 


( 2 . 11 ) 


For the present study, stationary grids (i.e., no body motion) are considered and 
the metric time derivative terms are zero. The transformation Jacobian and 
metrics have their geometric interpretations: the transformation Jacobian is the 
inverse of the local grid cell volume, and the metrics are grid cell area 
projections. They give good indication on the grid quality, for example, whether 
the grid is stretched too rapidly. 

The transformed Navier-Stokes equations can be written in strong 
conservation-law form (see, for example, Anderson et al [1984]) in a generalized 
coordinate system ( ^ , r], £, t) as follows: 
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where the symbol A indicates transformed variables. The transformed, 
conservative flow variables vector and the inviscid flux vectors are 


"p ' 


pu 

pu 


puU + £ x p 

pv 

, E = J~ l 

pvU + £ y p 

pw 


pwU + % z p 

e 


(e+p)U-Z tP 


pv 


pw 

puV + TJ xP 


puW+^p 

pvV + 77 y p 

, G = J~ l 

pvw+( r p 

pwV + Tj zP 


pwW + 

(e + p)V -T] t p_ 


(e+p)W-( lP 


( 2 . 13 ) 


( 2 . 14 ) 


where the variables U , V and W are the contravariant velocity components. The 
contravanant velocity U is the component of velocity parallel to the wing surface 
and in the direction of the wing chord, V is the component of velocity in the 
spanwise direction, and W is normal to the wing surface. They are related to the 
velocity components, u, v, and w, as follows: 

u = & + & + G,v+G z w 

V=T}' + ri x u+7i y v + Ti z w ( 2 . 15 ) 

*-C, + 0 + tv + $> 


The transformed viscous flux vectors are 


CHAPTER 2. GOVERNING EQUATIONS AND NUMERICAL METHOD 


20 


0 


E.. = J 


-i 


SA + ZyPy + ZA 


F=J 


-1 


G.. = J 


-1 


0 

*7,^ + 1?,^ + 17^ 
W yx + 7} y r„ + Tl z T yz 
Wv + Wv + W* 

. nA + nfiy + nA , 
o 

CxV + ^ T xy + C^z 

CA+CA+CA 


(2.16) 


(2.17) 


(2.18) 


2.1.3 Thin-Layer Approximation 

Most aerodynamic applications of practical interest, including the trailing vortex 
problem, involve high Reynolds number flows. For the high Reynolds number 
flows, the viscous effects are confined to the thin shear layers near the body 
surface, in the wake and in the trailing vortex core. Due to the limitations in 
computer memory and computational time it is necessary to concentrate the 
available grid points in these thin shear layers. To resolve the corresponding 
derivative components in the viscous terms, fine clustering of grid points has to 
be made across these thin shear layers. This results in grid spacing that is fine, 
normal and near to the surface, and that is relatively coarse, tangential to the 
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surface. With this type of grid, even if the full Navier-Stokes equations were 
solved, the viscous terms possessing velocity gradients tangential to the body 
would not be resolved because of insufficient grid density along the surface. For 
most cases of interest (i.e. high Re flows), however, these terms are negligible 
anyway. Therefore, it is justifiable to eliminate from the calculation the viscous 
fluxes associated with the directions parallel to the surface, i.e., the £- and re- 
directions. This approximation is easily applied since the equations are already 
transformed into the body-fitted computational domain. 


The thin-layer approximation is motivated by the success of the boundary 
layer theory. All the assumptions, except one, made in the boundary layer 
theory are adopted in the thin-layer approximation. In the boundary layer 
theory, pressure is assumed to be constant across the boundary layer. However, 
the thin-layer approximation is less restrictive: it retains the normal moment um 
equation and allows pressure variation across the boundary layer. 


Applying the thin-layer approximation (i.e., retaining only the d /d£ 
terms and eliminating all the d /<?£ and d /9 tj terms in the viscous stress 
components) then, the non-dimensional, three-dimensional, unsteady Navier- 
Stokes equations in conservation-law form in transformed, body-fitted, 
curvilinear coordinates become: 


d0 + 5£ + aF + <9G = 1 dS 
dr + + dr] + d£ ~ Re d£ 


(2.19) 


where 
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0 


s = r l 


/dm 1 Ur + (li/3)m 2 {; x 

mjn 3 + (^/3)m 2 ( 77 x m + 7 j y v + tj z w) 


( 2 . 20 ) 


with 


”h = $>c + £>f + 4>f 


( 2 . 21 ) 


m. 


= — (u 2 + v 2 + W> 2 ) H j— r(a 2 ) 

2 V ^ /V(y-l) v 


2.2 Numerical Algorithm 

An implicit, noniterative, time-accurate, finite difference scheme developed by 
Fujii and Obayashi [1986a, 1986b] is used to solve the three-dimensional 
compressible thin-layer Navier-Stokes equations. The scheme, which is first- 
order accurate in time and second-order accurate in space, employs the LU-ADI 
factorization. The basic structure of the algorithm is based on the Beam- 
Warming algorithm (see Beam and Warming [1976, 1978]) and incorporates the 
Steger- Warming flux vector splitting technique (see Steger and Warming [1981]) 
and diagonally dominant factorization (see Pulliam and Chaussee [1981] and 
Pulliam and Steger [1985]). The development of the algorithm is outlined below. 
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2.2.1 Beam and Warming Algorithm 

The Beam and Warming algorithm (see Beam and Warming [1976, 1978]) 
employs an implicit, approximately-factored, non-iterative technique. Implicit 
methods are preferred over explicit methods especially in obtaining the steady- 
state solution. A large time step can be used in implicit methods to reach the 
steady-state solution rapidly; whereas explicit methods suffer the disadvantage of 
having a severe restriction on time step size in order to maintain stability. 


A first-order accurate implicit time integration scheme is selected to 
march the solution of the unsteady Navier-Stokes equations in time. Second- and 
higher-order time accurate schemes are not used since they require the storage of 
the solution from previous time levels and result in a significant increase in the 
computer memory requirements. Discretizing Equation (2.19), using the first- 
order accurate backward Euler time integration scheme gives 


Q n+l -Q n +h 


dE n+l dF n+y dG n+l 
■ + — — + • 


1 dS 


>n+l 


dr\ Re d£ 


= 0 


( 2 . 22 ) 


A 

where h is the time step, n + 1 is the time level at which Q is desired, n is the 

A A A 

previous time level at which Q is known everywhere, and Q 1 = Q(nh). 


To have a non-iterative solution method, the flux vectors are linearized as 


follows: 
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£" +1 =£"+A" AQ'+0{h 2 ) 
F n+l =F n + B" AQ" + 0[h 2 ) 
G n+l =G” + C n AQ n + 0(h 2 ) 
S"* 1 =S” + M" AQ n +0(h 2 ) 


(2.23) 


where AQ* = Q H+l - O' and the flux Jacobian matrices A H , B”, C" and M n are 
given by: 






(2.24) 


These flux Jacobian matrices are derived analytically in Pulliam and Steger 
[1980]. 

The alternating direction implicit (ADI) algorithm was used by Beam and 
Wanning [1976, 1978] to replace the inversion of one huge matrix (which 
would be prohibitively expensive to compute) with the inversions of three block 
tridiagonal matrices, one for each direction, using a block tridiagonal solver. By 
applying the ADI algorithm. Equation (2.22) can be written in the “delta form” 
as follows: 

I + i b hS f A- + -KD,\,\ 

[ / + i, hS ( C" - i b hRe~ l 6 ( J-'M"J - i, D,| ; ] (2.25) 

= -i,h[s t ir + S,F ’ + Sfi" - +d c |„ + 

where I is the identity matrix. D, and D E are, respectively, the implicit and 
explicit artificial dissipation terms required for numerical stability. The “delta 
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form” is normally used for solving steady-state problems since the final solutions 

do not depend on the choice of implicit operators once steady state is reached 
(i.e., AQ" ->0). 

Notice that the solution blanking feature of the chimera scheme (see Benek 
et al [1985]) is incorporated into Equation (2.25) by including an integer i b 
which is assigned a value of zero or one at every grid point. If i b = 0, AQ" 
becomes zero and the solution remains unchanged at this grid point. If i b =1, the 
grid point is not blanked and is solved as part of the implicit solution. The 
usefulness of this blanking feature in modeling wing-tip blowing is discussed in 
Chapter 4. Fejtek and Roberts [1992] employed this feature for tilt-rotor 
computations. 

2.2.2 The Diagonally Dominant Factorization 

By applying the diagonally dominant factorization proposed by Pulliam and 
Chaussee [1981], the block tridiagonal matrices can be transformed into scalar 
tridiagonal matrices which are much cheaper computationally to invert. 
Recognizing that the flux Jacobian matrices A, B and C each has real eigenvalues 
and a complete set of eigenvectors, they can be diagonalized by similarity 
transformations as follows: 

A a =T-;'AT v A B = T;'BT n , A c = ^CT< (2.26) 

where A a , A b and A c are diagonal matrices containing the eigenvalues of 
matrices A, B and C, respectively. The elements of the diagonal matrices are the 
characteristic speeds of the flow. For example, 
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A. 


~U 0 

0 u 
0 0 
0 0 
0 0 


0 

0 

u 

0 

0 


0 

0 

0 

0 


0 

0 

0 

0 


U-aJg + l 



(2.27) 


Analytical expressions for the similarity transformation matrices 7^, T n and 7^, 

their inverse matrices, and the other diagonal matrices can be found in Pulliam 
and Chaussee [1981]. 


In the | -direction, for example, the Beam-Warming ADI operator can be 
written in the diagonal form as 


[/+i,A£ { A 




= T.Tl' + h hS f [T t K A T ?) - i„ T ( D , | { 77 ' 


(2.28) 


where the implicit smoothing factor e, = K,ha A with K, a user-specified 
constant, h the time step and o A the spectral radius of the matrix A. 

Moving 7^ and T~ l outside of the difference operator 8^ introduces an 

error which renders the method (at best) first-order accurate in time (see 
Pulliam and Chaussee [1981]). For steady-state calculations, the right hand side 
of Equation (2.25) goes to zero as A Q 1 — > 0, and the converged solution obtained 
using the diagonal algorithm is identical to that obtained from the original Beam- 
Warming ADI scheme since the right hand side is the same for both methods. 
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2.2.3 The LU-ADI Algorithm 

The scalar tridiagonal matrices of the implicit operators can be further reduced 
to products of lower and upper bidiagonal matrices (see Obayashi and Kuwahara 
[1986], and Fujii and Obayashi [1986a, 1986b]). This results in a further 
reduction in the computational cost. 


By decomposing the central differencing in Equation (2.28) into two one- 
sided differences using the flux vector splitting technique of Steger and Warming 
[1981], the £ -direction operator becomes 


I+i b hS f A-i b D t \ 4 


= T s [l+i b V f A\+i l A i A A ]r i ' (2.29) 


with 

A\^(A t ±\A,\)±r'e,J (2.30) 

where A\ contains all the positive eigenvalues and A\ contains all the negative 
eigenvalues, e, is the implicit artificial dissipation term, and J~ l is the inverse of 
the Jacobian evaluated at the central point 

Using first-order upwind differences. Equation (2.29) can be written as 

['H hS,A-i b A| ( ] = T ( [L a + D A + U A ]r f ‘ (2.31) 

where L A , D A and U A are lower bidiagonal, diagonal and upper bidiagonal 
matrices, respectively, and are given by 
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(2.32) 


Applying the diagonally dominant factorization (see Lombard et al [1983]) gives 
A. + D a + U A = (I, + Z>„ )D;‘(D X + U A ) + 0(A’ ) (2.33) 

where higher order terms are dropped since D A is of order 1 and L A and U A are 
of order h. 


Substituting Equation (2.33) into Equation (2.31), the LU factorization for 
the ADI operator in the £ -direction becomes 


/ + h hS ' A - 4 D, \ s ] = T s [ 4 + D a ] [p;‘(D,-H/,)] 77* 


lower bidiagonal upper bidiagonal 


(2.34) 


A similar procedure is followed for the tj- and ^-directions. The block 
tridiagonal matrix inversion for each direction has been reduced to a product of 
a lower and an upper scalar bidiagonal matrix. It is implemented by performing 
a forward sweep followed by a backward sweep. 


To ensure adequate stability of the thin-layer viscous terms while using the 
diagonally dominant factorization, it is required to add a small amount of 
additional dissipation to the split diagonal matrices A* (see Fujii and Obayashi 

[1986a]) as follows: 


A i c = ^(A c ±\A c \)±r , e,J±vI 


(2.35) 
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where 


Rep AC, 


(2.36) 


Finally, the LU-ADI algorithm can be summarized as follows: 

r s (l a + D t )d;'(d a + U A )r t % (l , + d.)d? (d,+u,)t;' 

T ; (Lc + D c )D?(D c + Vc)r ( 'AQ" 

= -i„ + s,r + S ; G ■ “*{$]- h [d £ | { + D E \ n + D e \ ( \q' 


(2.37) 


Analytical expressions for T~ l T n and T^T ; and their inverses can be found in 

Pulliam and Chaussee [1981] to reduce the computational effort. The inversion 
process in each direction consists of one forward scalar sweep and one backward 
scalar sweep. 


2.3 Turbulence Model 

For high Reynolds number flows, turbulence closure is required to account for 
the Reynolds stresses. By relating the Reynolds stresses to the rates of strain 
using the Bousinesq assumption, the effect of turbulence can be approximated by 
an eddy viscosity which accounts for the additional mixing caused by the 
turbulent flow. A two-layer algebraic turbulence model of Baldwin and Lomax 
[1978] is used in the present study. Since jet flows possess different turbulence 
characteristics, an eddy viscosity model proposed by Roberts [1987] for turbulent 
curved wall jet is implemented in the wall jet region for the wing-tip blowing 


cases. 
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2.3.1 The Baldwin-Lomax Model 


The eddy viscosity is evaluated in inner and outer layers, 

< 




M, 

(mA 


r ’ y — y crossover 
■ ’ y ^ crossover 


(2.38) 


where y is the normal distance from the wall and y aoaova is the smallest value of 
y at which the inner and outer values of fi t are equal. 


The eddy viscosity in the inner region is estimated using the Prandtl-Van 
Driest formulation: 

(dinner = M ( 2 ‘ 39 ) 

where |m| is the magnitude of the vorticity (|V x v|) and / is the “mixing” length 
scale given by 

/ = fcy [l - exp(- y + /A + )] (2.40) 

where k = 0.4, y + = y-Jp~ w T w /fi w with the subscript w denotes values at the wall 
and A + = 26. 

The Clauser formulation used in the Cebeci model for the outer layer is 
replaced by 

(2.41) 

where the Clauser constant K = 0.0168, C cp =1.6 and 

= min(y max , C^y^U 2 ^ / ) 


(2.42) 
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with C„* = 0.25 and U Sf is the difference between the maximum and the 
minimum total velocity magnitudes in the profile along the y-direction. The 
minimum total velocity magnitude is zero except in the wakes. The values of 
and are determined from the function 

F(y) = y|m|[l-exp(-y + M + )] (2.43) 

The quantity is the maximum value of F(y) in the profile along the y- 
direction, and y^ is the value of y at which it occurs. The function F Kleb (y ) is 
the Klebanoff intermittency factor given by 

^(3’) = [l + 5.5(C &t >/y_)‘f (2.44) 

with C KIeb = 0.3. 

The total effective viscosity can then be obtained as the sum of the laminar 
viscosity and the turbulent viscosity (n t ): 




(2.45) 


The laminar viscosity is determined from Sutherland’s formula: 


=Vr'f 





T^ + 198.6 °R 
T+ 198.6 ft 


where T is the temperature in degrees Rankine ( ft). 


(2.46) 


The total effective coefficient of thermal conductivity can be obtained by 


k = k t +k, 


SpEi | C JL 

Pr, Pr, 


(2.47) 
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For the range of temperatures and pressures of interest here, for air, the laminar 
Prandtl number Pr, =0.72 and the turbulent Prandtl number Pr i =0.90 (see, 

for example, Anderson et al [1984]). 


2.3.2 Turbulence Model for Wall Jet 


The algebraic eddy viscosity model developed by Roberts [1987] for turbulent 
curved wall jets is used in the wall jet region on the wing surface for the axial 
wing-tip blowing cases. This model has been previously used in the leading edge 
blowing for delta-wing (see Yeh et al [1989]) and tilt-rotor (see Fejtek and 
Roberts [1992]) applications. 


Assuming self-similar mean velocity profiles, the eddy viscosity of a wall 
jet is given by 


V,= 


JL b V 

4 k 2 m “ 




(2.48) 


where K = 0.073, k = 0.8814, is the maximum velocity, £ is the distance 

normal from the wall, b is the normal distance from the wall where the velocity 
is V /2, and C m , is the normal distance from the wall where the velocity is 
V . By correlating with the experimental results, b is found to have a value of 
about 7£ max . For £ > £/£«, is set to one. 


By including the effects of the wall curvature, the eddy viscosity is 
modified as follows: 


= — V C 

' max J max 


( 


\ b max J 


V/R ) 

dVM) 


(2.49) 
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where R is the radius of curvature of the wall. This wall jet model is applied 
from the jet exit slot to the flow separation point. 


2.4 Artificial Dissipation 


The Beam- Warming algorithm, although unconditionally stable for two- 
dimensional flows, is unconditionally unstable for three-dimensional flows (see 
Pulliam [1984]). Artificial dissipation is therefore necessary to stabilize the 
scheme. In addition to the explicit artificial dissipation, implicit artificial 
dissipation is added to increase the stability bound imposed by the explicit 
artificial dissipation and to speed up the convergence rate. 


The commonly used fourth-order artificial dissipation tends to produce an 
undesirable oscillatory solution near flow discontinuities such as shock waves. 
The second-order artificial dissipation, on the other hand, is too dissipative for 
other part of the flowfield even though it is able to damp out oscillations near 
flow discontinuities. An artificial dissipation model employing a nonlinear 
combination of second-order and fourth-order smoothing proposed by Obayashi 
et al [1988] is used in the present study. 


The explicit smoothing term, for the £ -direction, is given by 


d £ | ( 0- = v. 


(/-4> .)[ — 

' '**' 1 J Jm 


W-VA 


f 


*,jQr 


(2.50) 


'M 


where 
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<J> = matrix containing the flux limiter function 

^E ~ ^ E k 

K e = input constant 

o A = spectral radius of flux Jacobian matrix A 

=it/i +a V«. !+ € + £ 

The flux limiter function O varies from 0 to 1 depending the local flow 
gradient. For a relatively smooth solution, a value near one is used so that only 
the fourth-order dissipation is applied. Whereas for large flow gradients, a 
value near zero is used so that the second-order dissipation terms become 
dominant. More details on the flux limiter function are discussed in Obayashi et 
al [1988]. 




Chapter 3 


Grid Generation 


3.1 General Remarks 

Grid generation is an essential part of the overall computation of a fluid flow 
problem. Over the years various grid generation techniques have been 
developed and they can be broadly categorized into structured and unstructured 
grid generation methods. The structured grid generation methods can be further 
sub-divided into three main categories: (1) analytic methods, (2) algebraic 
methods, and (3) schemes based on partial differential equations. Each of these 
methods has its strength and weaknesses. These methods are discussed in more 
details in Anderson et al [1984]. 

To fully exploit the advantages of the grid generation method to be used, 
it is important to have a good understanding of the nature of the governing 
equations, the geometry of the flow field, and the expected flow behavior. For 
example, schemes based on partial differential equations are very attractive for 
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generating very smooth grids and allowing good grid spacing and orthogonality 
control. Their greatest drawback is that they are computationally expensive 
especially for three-dimensional applications. So far, most of the trailing vortex 
studies have employed geometrically simple lifting surfaces. Wings of 
symmetrical airfoil section, no twist, and no taper have been commonly used 
(see, for example, Spivey [1968], Yip and Shubert [1976], McAlister and 
Takahashi [1991], and Chow et al [1993]). The simple geometry allows a three- 
dimensional grid with good grid properties to be generated economically using 
schemes based on partial differential equations. This is achieved by first 
generating a two-dimensional grid at one spanwise station using schemes based 
on partial differential equations. A three-dimensional grid enveloping the 
complete lifting surface and the desired flow domain is then obtained by 
“stacking” the two-dimensional grid along the wing span. A grid generation 
method based on elliptic partial differential equations first proposed by 
Thompson et al [1974] is employed to generate the two-dimensional grid. Details 
of the method are presented in Section 3.2. 

Grid generation methods based on the unstructured grid philosophy have 
the advantage of tackling problems with complex geometry and flowfield 
elegantly. Strawn [1991] applied the unstructured adaptive-grid method to 
calculate the trailing vortex flowfield by solving the Euler equations. One 
difficulty encountered by Strawn [1991] was that the aspect ratio of the 
tetrahedral elements became more and more nonuniform for regions farther 
downstream. Flow solutions became deteriorated quickly in these regions. 
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3.2 Elliptic Grid Generation Method 

3.2.1 Thompson, Thames, and Mastin Method 

The ability of the elliptic grid generation method to generate very smooth and 
well behaved grids has made it one of the most widely used grid generation 
techniques for aerodynamic applications. The method was first proposed by 
Thompson, Thames, and Mastin [1974]. It was motivated by the smoothness 
properties associated with elliptic partial differential equation solutions. 
Successful applications by Thompson et al [1975, 1977a, 1977b] further 
popularized this method. 

In the present application, the elliptic grid generation method is used to 
generate a smooth, body-fitted, curvilinear two-dimensional grid. Grid points 
along the inner boundary, which includes the airfoil surface, and the outer 
boundary are specified. Interior grid points are first obtained by interpolation. 
They are subsequendy smoothed by the elliptic grid generation scheme which is 
oudined below. 

The mapping of the physical domain to the computational domain is 
defined by the Poisson equations as follows: 

$,+«==?(«) (31) 

where (jc,z) are the coordinates in the physical domain and (£,£) are those in 
the computational domain. The source terms, P and Q, are used to control the 
grid properties. Grid points in the computational domain, for the two- 
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dimensional application, are rectangular and evenly-spaced. Thus, in this way, 
(£, 0 are known. The unknowns are (x,z) in the physical domains. To 

determine (x,z) the dependent and independent variables in Equation (3.1) are 

interchanged. This results in the following transformed equations: 

(X x — 2(3 x + y X g ~ {p Q x^\ 

,, , (3.2) 

cc ~ 2(3 z ^ + y Zg — ~J yP z^+ Qz^ 


where 


a = x\ + z\ 

P = x i x l .+z i z i 
y=x\+z\ 


J = X;Zr - XfZf 


(3.3) 


Equation (3.2) can thus be discretized and solved numerically as long as P and Q 
are specified. For the present application, the alternating direction implicit 
(ADI) algorithm, which uses the approximate factorization to convert the left- 
hand-side matrix of Equation (3.2) into two tridiagonal matrices to enhance 
computational efficiency during matrix inversion, is employed. Details of the 
implementation of the algorithm can be found, for example, in Holst [1983]. 


3.2.2 Middlecoff-Thomas Algorithm 

It is noticed in Equation (3.1) that when P = Q = 0 the Poisson equations are 
reduced to the Laplace equations. This will provide no control over the grid 
point spacing near a boundary. The grid points tend to be pulled away from the 
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surface (see, for example, Sorenson and Steger [1980]). The specification of P 
and Q therefore becomes a main concern. The original P and Q terms proposed 
by Thompson et al [1974] require four user-specified coefficients to control the 
grid properties. The major setback of this proposal is that these user-specified 
coefficients are rather difficult to obtain in a relatively automatic way. 

To overcome this difficulty, Steger and Sorenson [1979], and Middlecoff 
and Thomas [1979] devised some automatic grid control techniques that 
significantly reduce the number of free parameters that need to be specified. In 
the present application, the Middlecoff-Thomas algorithm, which is outlined 
below, is employed. 

The P and Q terms are defined by Middlecoff and Thomas [1979] in 
terms of two other functions <f> and y / , which are given by 

*>=*(«)(£+ S’) 

G=r(£.4:)(£+C) 

Substituting Equation (3.4) into Equation (3.2) gives 

“(*« +0X')-2P x fl +r(x«+ r* { )=° 

“(z {! + )-2/3z, ; + r(z K + Vz ; )=0 

The 0(£,£) and terms are established in an automatic way from user- 

specified boundary conditions on x and z as follows: 


For £ = constant boundaries. 
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For £ = constant boundaries. 


0(£.f«.)= 



for 

kkk 


~ Z ttl 

u\ 

for j 

z fkk 

il 

~ X ttl 


for 

hkh 


~ z k! 

US 

for 

^kk 


( 3 . 6 ) 


( 3 . 7 ) 


All derivatives are discretized using central differences. Interior values of 
0 and iff are obtained using linear interpolation. 

The Middlecoff- Thomas algorithm has been proven to be simple and 
effective in grid control. Furthermore, interior grid point distribution tends to 
approximately mimic the boundary distribution of grid points. Grid points that 
are clustered in one direction on the boundary will tend to re main clustered in 
that same fashion in the grid interior. This allows the user a better control in 
clustering grid points to capture regions of high flow gradient. 
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3.3 Grid Details 

3.3.1 Geometry of Wing 

To illustrate the grid generation used in the present trailing vortex study, a semi- 
span rectangular wing with a constant and untwisted NACA 0015 symmetric 
airfoil profile along the entire span and a rounded wing-tip shape is used. The 
wing has a full-span aspect ratio of 6.75. The same geometry was used in the 
fixed-wing trailing vortex experiment of McAlister and Takahashi [1991], which 
was conducted in the NASA Ames 7- by 10-Foot Subsonic Wind Tunnel No. 2. 
A perspective view of the semi-span wing is shown in Figure 3.1. 

3.3.2 Two-Dimensional Grid 

Due to the simplicity of the geometry, a two-dimensional grid is first generated 
using the elliptic grid generation method (discussed in Section 3.2). Symmetry 
of the airfoil section also enables only the upper half of the grid to be generated 
(see Figure 3.2). This further reduces the cost in grid generation. 

As discussed in Section 2.1, the application of the thin-layer 
approximation to the governing Navier-Stokes equations has limited the 
computer memory and computational time requirements to a manageable level. 
To resolve the viscous fluxes normal to the lifting surface, which are the 
dominant contribution to the overall viscous fluxes, very fine grid points must be 
clustered near the lifting surface so as to capture the high velocity gradient in the 
boundary layer. In addition, very fine grid points must be clustered aft of the 
trailing edge in order to capture the wake, the lifting vortex sheet, and the 
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trailing vortex effects. The C-grid topology, as opposed to the O- and H-grid 
topologies, is selected because grid points can be clustered more efficiently to 
capture these flow features. Furthermore, the O-grid topology has a mapping 
singularity at the trailing edge of all sharp-trailing-edge airfoils; and the H-grid 
topology places a singularity at the leading edge of all blunt-leading-edge 
airfoils. The C-grid topology, on the other hand, is free from these singularity 
mapping problems. 

Fine grid clustering in the streamwise direction is also required at the 
leading edge of the airfoil section to capture the high velocity gradient and the 
high leading edge suction peak. The trailing vortex problem also demands fine 
grid clustering at the trailing edge in the streamwise direction, especially near 
the wing-tip region. It was observed (see, for example, Spivey [1968], Yip and 
Shubert [1976], McAlister and Takahashi [1991], and Chow et al [1993]) that 
high suction peak (or peaks) occurred on the upper surface of the wing-tip 
region due to the presence of the rolled-up tip vortex (or vortices depending on 
the geometry of the wing-tip). To capture this effect, fine grid has to be 
clustered there. 

Aft of the trailing edge near the wing tip region, the trailing vortex 
continues to develop and roll up. The distance required for the trailing vortex to 
completely roll up depends on the aspect ratio, lift distribution, etc, of the wing 
(see Spreiter and Sacks [1951]). Fine grids must therefore be clustered aft of the 
trailing edge near the wing tip region for some downstream distance to capture 
the formation of the trailing vortex. 
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To ensure fine and smooth grid clustering in the streamwise direction near 
the leading and trailing edges of the airfoil section, a cosine stretching function is 
employed. Stretching is also employed aft of the trailing edge in the downstream 
direction to ensure fine clustering near the trailing edge and smooth stretching of 
grid farther downstream. It is found to be difficult to use the cosine stretching 
for grid points aft of the trailing edge. Smooth grid transition near the trailing 
edge cannot be obtained with the present amount of grid points available because 
the downstream distance is much greater as compared to the chord length. A 
geometric stretching is used instead to cluster grid points aft of the trailing edge. 

The geometric stretching is probably the simplest stretching method. To 
distribute N number of points along a curve of length 5, with the arc length 
between the first two points of AS , the total length can be written as 

AS + oc AS + cc 2 AS 4- ... + cc N 2 AS + ct N *AS — S (3.8) 

where a is the stretching factor. 

The Newton-Raphson iterative root finding method can be used to 
determine a such that Equation (3.8) is satisfied. It is implemented by 
multiplying Equation (3.8) by a and then subtract this new equation from 
Equation (3.8). This results in the following equation: 

(a-l)S-(a*-l)AS = 0 = / (3.9) 

Notice that a function / is defined in the iterative procedure. The objective of 
the iteration is to determine a value a such that / = 0 is satisfied within a 
desired tolerance, a is determined iteratively via 
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«„ + i = a„— (3.10) 

where n is the iterative time step and 

f' = S-Na N ~ l AS (3.11) 

This method does not provide any control on the size of the last grid spacing. A 
more sophisticated method developed by Vinokur [1983] allows control on the 
last grid spacing. However, for external flow aerodynamic applications, the 
control of the size of the last grid spacing is not critical. Geometric stretching is 
also used to generate grid points normal to the wing surface. 

3.3.3 Three-Dimensional Grid 

To generate a three-dimensional grid, the two-dimensional grid, illustrated in 
Figure 3.2, is “stacked” along the wing span starting from the upper surface of 
the wing center-plane. The two-dimensional grid is then wrapped around the 
wing-tip to form a rounded wing-tip. “Stacking” of the two-dimensional grid 
continues until the lower surface of the wing center-plane is reached. 
Orthogonality of the two-dimensional planes with the wing surface is ensured 
during the “stacking” process. A cut-away- view in the y-z plane of the grid 
generated is shown in Figure 3.3. (The standard right-hand aeronautical 
convention is used for the x-y-z body-axis coordinate system: the origin of the 
body-axis coordinate system is at the intersection of the wing center-plane and 
the wing leading edge, x is positive from the leading edge to the trailing edge, y 
is positive toward the tip of the pilot’s right wing, and z is positive upward.) 
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Similar to the experimental setup, only the semi-span of the wing is 
considered due to the symmetry of the flow configuration. Implicit in the use of 
a plane of symmetry is the assumption that the flowfields on either side of the 
center-plane are a mirror image of each other. To facilitate the specification of 
the boundary condition at the center-plane, the two-dimensional grid actually 
extends one grid plane inboard of the plane of symmetry; and the distances of 
these two planes from the plane of symmetry are equal. Flow conditions at this 
plane are forced to be those at the first plane from the center-plane on the pilot’s 
right wing to preserve symmetry of the flow. This arrangement is deemed to 
simulate the actual flow more closely as compared with the experimental setup 
whereby boundary layer growth from the wall at which the wing is mounted 
may interfere with the actual flow. 

Cosine stretching is used to cluster grid points in the spanwise direction. 
This allows fine clustering of grid points around the wing tip so as to capture the 
trailing vortex better. A cut-away isometric view of the wing and the three- 
dimensional computational grid is shown in Figure 3.4. 

Most of the results presented are based on computations using 83x69x61 
(i.e., 349,347) grid points. There are 83 points in the | -direction with 48 points 
defining half of the airfoil section and 35 points downstream of the trailing edge. 
The grid extends 8 wing chords from the trailing edge to the outflow boundary. 
69 two-dimensional grid planes are “stacked” in the spanwise T ] -direction with 
25 planes defining the rounded wing tip. There are 60 grid points stretch from 
the wing surface to the outer boundary in the ^-direction. The grid extends to 6 
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wing chords normal to the wing surface and 6 wing chords ahead of the wing 
leading edge. 
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Figure 3.1: Geometry of the semi-span wing having a NACA 0015 symmetric 
airfoil section, with no twist, no taper, no sweep, and a rounded wing-tip. 
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Figure 3.2: Two-dimensional grid at the wing center-plane showing the C-grid 
topology. Fine clustering of grid points at the leading edge, the trailing edge, aft 
of the trailing edge, and close to the wing surface. 
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Figure 3.3: A y-z plane view (looking upstream) of the three-dimensional grid 
showing the rounded wing-tip of the wing and the fine clustering of grid points 
around the wing-tip to capture the wing-tip vortex. 
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Figure 3.4: A cut-away isometric view of the wing and the three-dimensional 
computational grid. 


Chapter 4 


Boundary Conditions 


4.1 General Remarks 

Boundary conditions at the flowfield domain boundaries are the driving force in 
determining the overall flowfield solution. In the implementation of the 
numerical algorithm, the boundary conditions can be specified either implicitly 
or explicitly. When specified implicitly, the boundary conditions are coupled 
and coded into the numerical algorithm. Explicit specification of the boundary 
conditions, on the other hand, requires that flow variables at the boundaries to be 
evaluated using the most recent solution. Due to the fact that explicit treatment 
of the boundaries leads to a far more simple and flexible scheme, where 
boundary conditions become a modular element that can be put in or pulled out 
of a computer program without disturbing the implicit algorithm (see Pulliam 
and Steger [1980]), it is implemented in the present study. The flexibility of the 
explicitly specified boundary conditions also allows wing-tip blowing (to be 
discussed in section 4.3) to be incorporated into the numerical study easily. 
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4.2 Boundary Conditions 

On the surface of the wing, no-slip boundary condition is imposed, that is, all 
components of velocity are set to zero (« = v = w = 0) for viscous (Navier- 
Stokes) computations. 

The pressure on the wing surface is obtained by solving the normal 
momentum equation (see Steger [1977] and P ulliam and Steger [1980]): 

p ( (W, + + S,f, )■ + p„{nA, + n,i, + n«?,)+p c (£ + C + C) 

= ~PV( (,u ( + i,v f + f > { ) - pV{i,u n + f,v, + ) (4 , } 

= 0 

Notice that p n = dp/dn = 0 at the body surface and n is the direction normal to 
the body surface. 

By setting U and V to zero at the body surface for viscous computations 
and discretizing the derivatives using second-order central differences in the £ - 
and T] -directions and second-order one-sided differences in the ^-direction, p at 
the body surface can be calculated by solving the implicit equations using the 
approximate factorization technique. 

Once the pressure at the body surface is determined, the total energy per 
unit volume, e, can be calculated using Equation (2.7). By ass umin g adiabatic 
wall conditions (i.e., no heat flux across the wall or dTjdn = 0), and using the 
relation dp/ dn = 0 at the wall, it is necessary that dp/dn = 0 at the body surface 
as well from the ideal gas equation (p = pRT). The density at the body surface is 
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then calculated using the zero-order extrapolation from the nearest point normal 
to the surface. 

Due to the symmetry of the flowfield only half of the wing span is 
modeled to reduce the computational cost. The symmetric boundary conditions 
at k = 1 and k = 3 (which are next to the plane of symmetry or mid-span: k = 2) 
are imposed as follows: 

A=Pa 

(pu\ = (pu) 3 

(pv) 1 =-(pv) 3 (4.2) 

(pw\ = (pw\ 
e x — £ 3 

For the outer boundaries, freestream conditions are imposed on the inflow 
boundary whereas characteristics-based formulation is applied at the outflow 
boundary. For example, the present trailing vortex problem involves subsonic 
flow in three dimensions. In this case, only one state variable is to be specified at 
the outflow boundary and the remaining four variables must be obtained from 
the characteristic relations because four of the characteristic velocities are 
positive and the fifth one is negative. Zero-order extrapolation is used in 
obtaining the four variables at the outflow boundary from the interior solution. 

4.3 Wing-Tip Blowing 

The blowing intensity is quantified by the blowing momentum coefficient, C li , 


which is defined as 
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ttpi«Vj ja -dA ja 



(4.3) 


where p ja is the jet exit density, V ja is the jet exit velocity, A Jel is the jet slot 
area, q_ is the freestream dynamic pressure, and is the wing reference area, 
which is simply the full-span wing planform area. 


43.1 Actuator Plane Concept 

For the case of axial wing-tip blowing, a step in the wing surface is formed by 
the blowing slot. To simply the grid generation required in this study, the jet 
slot is not resolved and an actuator plane concept is employed to model the wall 
jet. A schematic sketch of the actuator plane concept is presented in Figure 4.1. 
The actuator plane is, in effect, a discontinuity imposed at the jet slot location 
whereby the flow variables undergo a discontinuous change and the jet is 
modeled as a one-sided source of mass, momentum and energy. This concept has 
been successfully applied to the F- 18 forebody (see Tavella et al [1990]), an 
ogive cylinder (see Font and Tavella [1991]), a tilt-rotor (see Fejtek and Roberts 
[1992]) and a delta-wing (see Craig [1993]). Tavella et al [1990] compared the 
results obtained by using the actuator plane concept with those obtained by Yeh 
et al [1989] who resolved the jet slot and found that the differences were not 
significant. 

For the case of spanwise wing-tip blowing, no step on the wing surface is 
formed by the jet slot. The jet can therefore be modeled at the wing surface in a 
straight-forward manner. 
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4.3.2 Solution Blanking Feature of Chimera Scheme 

In modeling the wall jet for the axial wing-tip blowing case, L je[ number of grid 
points normal to the wing surface at the jet slot location are used to define the jet 
slot thickness. In order to avoid the flowfield solution being solved at the jet slot 
location, Craig [1993] used the multi-zonal approach and defined the jet slot 
location at the inter-zonal boundary. This approach is not very practical because 
when one changes the jet slot location, for example for parametric study, one has 
to change the entire computational grid and re-calculate the entire flowfield. 
The approach of employing the solution blanking feature of the chimera scheme 
adopted by Fejtek and Roberts [1992] is instead used in the present study. By 
blanking off the implicit solution at the jet slot location and updating the solution 
explicitly with the “interior boundary conditions” defined by the jet, the jet slot 
can be located anywhere in the flowfield. The allows great flexibility in locating 
the jet slot without confining it to the inter-zonal boundary. 

The static pressure at the jet slot exit is assumed to be the local static 
pressure just outside the jet width and constant across it for the axial wing-tip 
blowing case. For the spanwise wing-tip blowing case, the static pressure at the 
jet slot exit is also assumed to be the local static pressure; but it is allowed to 
vary across the jet width. Wing-tip blowing is applied by changing the total jet 
plenum supply pressure which is normalized by the freestream ambient pressure, 

i* e -’ P plenum /P~ ‘ 

The temperature at the jet slot exit is assumed to be freestream ambient. 
The jet exit density, p jet , is calculated using the ideal gas equation. By assuming 
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isentropic expansion of the compressed air from the plenum pressure, p pUmon , to 
the local static pressure, p, the jet exit velocity is obtained by 




\-(r-0/r 


PpUmon 


-1 


2 r p 


r-ip 


(4.4) 


The total energy per unit volume at the jet slot exit is determined by Equation 
(2.7). 








Chapter 5 


Results and Discussion 


5.1 General Remarks 

Results obtained in this study are presented and discussed in this chapter. The 
chapter is divided into no blowing cases, axial wing-tip blowing cases and 
spanwise wing-tip blowing cases. Computational results are compared with 
experimental data and theoretical analyses. The effects of axial and spanwise 
wing-tip blowing on the performance of the wing and the behavior of the trailing 
vortices are also discussed. 

5.2 No Blowing Cases 

5.2.1 Configuration 

The configuration used in the present study for no blowing and axial wing-tip 
blowing cases is the fixed-wing trailing vortex experiment of McAlister and 
Takahashi [1991]. The experiment was conducted in the NASA Ames 7- by 10- 
Foot Subsonic Wind Tunnel No. 2. Semi-span configurations were considered in 
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the study and all the wings were rectangular and had a constant and untwisted 
NACA 0015 symmetric airfoil profile along the entire span. The lateral edge of 
each wing-tip was machined to a flat or square face, but could be made round by 
the addition of an end cap. For the present study, the rounded wing-tip 
configuration is considered. A perspective view of the semi-span wing is shown 
in Figure 3.1. With the rounded wing- tip end cap, the aspect ratio of the wing is 
6.75. Three different chord lengths of 12.0, 16.2 and 20.4 inches were used in 
the experimental study. However, the aspect ratio remained fixed by changing 
the length of the wing span. 

Computations were performed for angle-of-attack a = 4°, 8°, and 12°, 
ffeestream Mach number M m =0.17, and freestream Reynolds number based on 

the wing chord Re = 2.0 x 10 6 . 

5.2.2 General Flow Features 

The rolling-up of the tip vortex and its subsequent convection downstream is 
depicted by the particle traces in Figure 5.1. It is observed that farther inboard, 
the flow remains very much two-dimensional. 

Figure 5.2 (a) shows the side view of the trailing vortex. It is noticed that 
the tip vortex starts to roll up even before it reaches the trailing edge of the wing 
and that it eventually rolls up into the freestream direction and is convected 
downstream. This is also observed in the experiments (see, for example, 
McAlister and Takahashi [1991] and Chow et al [1993]). Figure 5.2 (b) shows 
the plan view of the trailing vortex. It is observed that as the trailing vortex 
moves downstream, it moves inboard asymptotically as well. For the present 
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rectangular wing, the trailing vortices move about 5% inboard at stations far 
downstream of the wing. Spreiter and Sacks [1951] showed that the trailing 
vortices could move inboard by as much as about 25% at stations several chord 
lengths behind the wing for an elliptically loaded wing. For elliptically loaded 
wings, the asymptotic spacing between the trailing vortices is independent of the 
angle of attack; whereas, for non-elliptically loaded wings, the asymptotic 
spacing is a function of the angle of attack (see Spreiter and Sacks [1951]). 

Figure 5.3 shows the velocity vectors at the quarter-chord location. 
Notice that the flow accelerates from the higher pressure lower surface towards 
the lower pressure upper surface by going around the wing-tip. No-slip 
boundary conditions at the wall are observed for the Navier-Stokes 
computations. 

5.2.3 Pressure Distribution 

The computational results for pressure (in terms of the coefficient of pressure, 
C p ) on the wing surface for the outer 25% of the wing span and for cases of 

a = 4°, 8°, and 12° are given in Figures 5.4, 5.5, and 5.6, respectively. The 
corresponding experimental results obtained by McAlister and Takahashi [1991] 
are also included in the figures for comparison purposes. The freestream flow 
conditions are: A/. = 0.17 and Re = 2.0 x 10 6 . 

As discussed in Section 5.2.1, the experiment was conducted in the NASA 
Ames 7- by 10-Foot Subsonic Wind Tunnel No. 2. Due to the presence of the 
wind tunnel walls, pressure measurements are altered from their free-air values 
because of blockage and distortion of the streamlines. In a closed test section, 



CHAPTER 5. RESULTS AND DISCUSSION 


61 


blockage has the effect of producing a more dense flow and a higher velocity in 
the region where the wing is located (see McAlister and Takahashi [1991]). The 
nondimensional pressure coefficients are based on the “freestream” static and 
dynamic pressures that are obtained from a pitot-static probe placed upstream in 
the test section. By taking into account the blockage effects, the corrected 
pressure coefficient is given as 



P~(P-* + AP-) 

<l~u+Aq_ 


(5.1) 


where the subscript u denotes an uncorrected value and the symbol A stands for 
the difference between corrected and uncorrected values. The corrected 
freestream velocity, which is more representative of the constricted flow in the 
test section where the wing is actually positioned, is given by 

V>(l + e)V„ (5.2) 

where i s the “freestream” velocity at an upstream location that is not 
influenced by the wing (measured with an upstream pitot-static probe) and e 
represents the increase of velocity due to blockage effects. 

By assuming constant density of flow due to the low Mach number and 
using the approximation (1 + e)~ 2 ~ (1 - 2e) , Equation (5.1) can be simplified as 

C p = C pu (l-2e)+2e (5.3) 

where is the pressure coefficient that would be formed using upstream 
reference values ( p _ and <?_), without regard for blockage effects. 
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Since computations were carried out based on “free-air” formulations, the 
wall correction procedures proposed by McAlister and Takahashi as described 
above are applied to the experimental data for consistent comparisons. It is 
observed in Figures 5.4, 5.5, and 5.6 that the computational results compare very 
well with the experimental data except at the region near the leading edge and at 
the trailing edge region near the wing-tip. The computations did not capture the 
sharp leading edge suction peak and the vortex induced low pressure peak near 
the trailing edge at the wing-tip region too well, especially for the higher angle- 
of-attack case, due to the computer memory constraint. Results can be improved 
by clustering more grid points in these regions. 

It is also noticed that near the wing-tip, tip vortex is formed on the upper 
surface of the wing. This tip vortex, which is reported by McAlister and 
Takahashi [1991] as well as Spivey [1968], Yip and Shubert [1976], and Chow et 
al [1993], are captured quite well in the computations. A single vortex induced 
low pressure peak is observed near the trailing edge of the upper surface for the 
rounded wing-tip configuration. Multiple pressure peaks have been reported for 
squared wing-tip configurations (see, for example, Spivey [1968], and McAlister 
and Takahashi [1991]). The formation of the tip vortex is clearly shown in 
Figure 5.7 in terms of the pressure contours and Figure 5.8 in terms of the 
velocity vectors. Low pressure is observed at the vortex core due to the high 
rotational velocity. Figure 5.8 also shows the highly three dimensional nature of 
the flow near the wing-tip. 

Further inboard of the wing span, the flow remains very much two- 
dimensional for most of the wing span. As illustrated in Figure 5.9, the C p 
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distributions for the inboard 60% of the wing span are very s imil ar This is also 
observed by McAlister and Takahashi [1991] though no experimental data are 
available for comparison. 

Figure 5.10 shows the stagnation pressure coefficient contours across the 
trailing vortex for various axial stations downstream of the tr aili ng edge. Notice 
that the trailing edge is at x/c = 1. The low stagnation pressure vortex core is 
clearly seen in the figure. It is observed that as the trailing vortex moves 
downstream, the vortex core size increases. 

5.2.4 Forces and Moment 

The computed coefficients of lift, drag and pitching moment (about the wing 
quarter-chord) for a = 4°, 8°, and 12°, Af_= 0.17 and Re = 2.0 x 10 6 are 

compared with the experimental data of McAlister and Takahashi [1991] in 
Figure 5.11. 

As discussed in Section 5.2.3 and in McAlister and Takahashi [1991], the 
presence of the wind-tunnel walls changes the flowfield around the wing. A 
change in streamline curvature (caused by the airfoil images due to the presence 
of the wind-tunnel walls) has the effect of imparting greater “apparent” camber 
to the airfoil and inducing a higher angle of attack (or an increase in the effective 
airfoil incidence). In order to have consistent comparisons with the computed 
results, the lift coefficients obtained in the experiments were corrected for wall 
effects as recommended by McAlister and Takahashi [1991]. For the present 
experimental setup, a 0.51° correction in the angle of attack is required for the 
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case of a = 12°. Effects of the wall on the coefficients of drag and pitching 
moments are deemed to be small; thus, no wall corrections are made. 

Overall, the computed coefficients of lift and pitching moment compare 
favorably with the experimental results. Drag coefficients, however, are over 
predicted by the computations. 

5.2.5 Vortex Strength 

Figure 5.12 shows the streamwise distribution of the normalized total circulation 
for the case of oc = 12°, Af_ =0.17 and Re = 2.0 x 10 6 . The vortex strength or 

the total circulation is defined as 

r = jv ds (5.4) 

s 

where V is the velocity vector, ds is the vector element of length along the path 
of integration, and s is the path along which the line integral is taken. Notice that 
r is normalized by cV_/2, where c is the wing chord and V. is the freestream 

velocity, so that direct comparison with the section lift coefficient at the mid- 
span of the wing can be made. 

According to the inviscid analysis, the strength (at either subsonic or 
supersonic speeds) of one of the trailing rolled-up vortices at stations far behind 
the wing is equal to the sum of the strengths of all the vortices shed from one- 
half of the wing, and, hence, it is equal to the magnitude of the circulation 
around the wing in the plane of symmetry (see Spreiter and Sacks [1951]). Thus, 
the normalized total circulation at stations far behind the wing is equal to the 
section lift coefficient in the plane of symmetry. 
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For stations immediately downstream of the trailing edge of the wing, the 
trailing vortices are not fully rolled-up and much of the vortex strength of the 
vortex system resides in the vortex sheet As the vortex system moves farther 
downstream, the vortex sheet finally rolls up completely into the trailing vortices 
(see Spreiter and Sacks [1951]). 

The computed normalized total circulation is compared with the result of 
the inviscid analysis in Figure 5.12. Notice that the path of integration for the 
circulation calculation was chosen to enclose the trailing vortex but exclude the 
vortex sheet as much as possible so that the rolling-up process of the vortex sheet 
can be captured. It is observed in Figure 5.12 that the normalized total 
circulation reaches a peak value at some station aft of the trailing edge (at which 
the vortex sheet has completely rolled-up) and then reduces its strength slightly 
until an asymptotic value is reached farther downstream. This asymptotic value 
is found to be about 85% of the theoretical value predicted in the inviscid 
analysis. The reduction in the circulation strength could be due to viscous 
dissipation as was also observed in the experimental studies of Higuchi et al 
[1986]. 

The comparison of the computed normalized downwash distribution 
downstream of the trailing edge with the theoretical analysis of Spreiter and 
Sacks [1951] for the case of a = 12°, = 0.17 and Re = 2.0 x 10 6 is shown in 

Figure 5.13. Notice that the downwash w is normalized by V^C L b/nARb' 
(where C L is the lift coefficient of the wing, b is the full-span of the wing, AR is 

the aspect ratio of the wing, and b' is the distance between the trailing vortices ) 
as used by Spreiter and Sacks [1951]. Since the distance between the trailing 
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vortices decreases as they move downstream, the spanwise distance, y, from the 
plane of symmetry of the wing at various downstream stations is normalized by 
b'/2. 

The inviscid analysis of Spreiter and Sacks [1951] assumed that the core 
size of the trailing vortices remained unchanged (i.e., the vortex core did not 
grow in size) and that the circulation strength remained constant; hence, a 
universal curve for the normalized downwash distribution for stations far 
downstream from the wing was obtained. In the present viscous computations, it 
is observed that although the circulation strength remains constant at various 
streamwise stations farther downstream from the wing, the maximum magnitude 
of the downwash (or upwash depending on whether it is inboard or outboard of 
the trailing vortices) decreases as the trailing vortices travel downstream because 
the size of the vortex core increases due to the entrainment of the surrounding 
fluid into the viscous vortex core. 

5.3 Axial Wing-Tip Blowing Cases 

5.3.1 Configuration 

The same configuration described in Section 5.2.1 for the no blowing cases is 
used for the axial wing-tip blowing cases. In addition, a jet slot is created as 
shown in Figure 5.14. The jet slot has a length of 11.2% of the wing chord and 
a width of 1.5% of the wing chord. Due to the presence of the jet slot, a small 
step of height 1.5% of the wing chord is created on the wing surface. As 
discussed in Section 4.3 the actuator plane concept is employed in the 
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computations to model this jet slot. The intensity of the axial wing-tip blowing is 
controlled by changing the jet plenum supply pressure. 

Results presented here are for a = 8°, =0.17 and Re = 2 x 10 6 . The 

blowing intensity varies from =0.0000 to 0.0047. The highest blowing 

intensity gives a jet exit velocity closed to sonic condition. Results presented are 
also for jet slot located at the 90%-chord of the wing. Computations were 
performed for jet slot located farther upstream as well; however, it was observed 
that less penetration into the flowfield was achieved as the wall jet had to 
overcome additional wing surface area. They are therefore not presented. 

5.3.2 General Flow Features 

Figure 5.15 shows the closed- up views of the velocity vectors in the wing-tip 
region for cases with the axial wing-tip blowing off and on. Without blowing, a 
typical turbulent boundary layer velocity profile is observed in Figure 5.15 (a). 
By turning the blowing on, the presence of the jet transforms the boundary layer 
profile into a wall jet profile (see Figure 5.15 (b)). It is observed that consistent 
with experimental observations the velocity profiles upstream of the jet slot are 
not affected by the presence of the jet. This consistency is possible since the 
actuator plane concept models the jet as a one-sided source of mass, momentum 
and energy. 

Figure 5.16 shows the effect of axial wing-tip blowing on the flowfield 
downstream of the trailing edge of the wing. As shown in Figure 5.16 (b), the 
high velocity jet penetrates into the vortex core of the trailing vortex for some 
downstream distance before the effect of the blowing dimini shes It is also 
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observed that due to the rolling-up process of the vortex sheet, the jet re-aligns 
itself into the freestream direction together with the trailing vortex. 

5.3.3 Vortex Strength 

Figure 5.17 shows the variation of the streamwise distribution of the normalized 
total circulation with the blowing momentum coefficient. As discussed in Section 
5.2.5 for the case with blowing off, the normalized total circulation reaches a 
peak value at some station aft of the trailing edge (at which the vortex sheet has 
completely rolled-up) and then reduces its strength slightly until an asymptotic 
value is reached farther downstream. With blowing on, the same phenomenon is 
observed except that now it takes longer for the circulation to reach a peak value, 
i.e., it takes longer for the vortex sheet to roll up completely. As the blowing 
intensity increases, the rolling-up process is delayed further. The delay in the 
rolling-up process due to blowing can be seen more clearly by the particle traces 
in Figure 5.18. It is shown in Figure 5.18 (a) that without blowing the tip vortex 
starts to roll up even before it reaches the trailing edge of the wing. With 
blowing turned on in Figure 5.18 (b), the momentum of the jet delays the 
rolling-up of the vortex sheet. 

It is also observed in Figure 5.17 that the total circulation strength of the 
trailing vortex increases slightly for stations farther downstream of the wing as 
the blowing intensity increases. This is due to the fact that the axial wing-tip 
blowing increases the lift on the wing slightly (to be discussed in Section 5.3.4), 
and that the higher the wing lift, the stronger the circulation of the trailing 


vortex. 
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Figures 5.19 and 5.20 show the downwash distribution at various axial 
stations downstream of the trailing edge for cases with blowing off and on, 
respectively. For stations immediately downstream of the trailing edge, the 
downwash differs slightly for cases with blowing off and on. However, the 
downwash distribution for the case with blowing on is quite similar to that with 
blowing off for stations farther downstream of the wing except that the former 
has a slightly higher value due to the stronger circulation strength discussed 
earlier. 

Since the circulation strength of the trailing vortex is an overall feature of 
the wing loading, the axial wing-tip blowing, although possibly able to change 
the behavior of the local flowfield, is unable to alter the flowfield features 
farther downstream. It was also observed in the experimental study of Snedeker 
[1972] that the axial wing- tip blowing, although it altered the near-field velocity 
distribution, did not reduce the rolling moment imposed on an object farther 
downstream. Dunham [1976] classified the axial wing-tip blowing as an 
unsuccessful concept for aircraft wake vortex minimiz ation. 

5.3.4 Forces and Moment 

Figure 5.21 shows the effects of axial wing-tip blowing on the overall 
performance of the wing in terms of lift, drag and pitching moment (about the 
wing quarter-chord). It is observed that as the blowing intensity increases the 
lift on the wing increases. This is due to the fact that the high velocity jet 
entrains the surrounding fluid on the upper surface of the wing. This creates a 
low pressure region in the vicinity of the jet (see Figure 5.22) and increases the 
lift on the wing. 
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Since the jet slot is located near the trailing edge of the wing, this low 
pressure region also increases the pressure drag on the wing and results in higher 
overall drag. Similarly, the low pressure region also creates a pitching down 
moment on the wing. Thus, the overall pitching moment (taken as positive for 
pitching up moment) decreases as the blowing intensity increases. 

5.4 Spanwise Wing-Tip Blowing Cases 

5.4.1 Configuration 

The configuration used in the experiments of Tavella et al [1985, 1988] and Lee 
et al [1989] (see Figure 5.23) is used in the present computational study with 
spanwise wing-tip blowing. The experiments were conducted in the Stanford 
low-speed wind-tunnel which had a 45.7 cm x 45.7 cm test section and a 
freestream velocity of 40 m/s, giving a chord-based Reynolds number of 
4 x 10 s . A semi-span rectangular wing with a constant and untwisted NACA 
0018 symmetric airfoil profile along the entire span was used. The wing-tip is 
rounded and the (full-span) aspect ratio of the wing is 3.28. 

A jet slot with a thickness of 0.16 cm was positioned in the plane of 
symmetry of the wing-tip, extended over 73.3% of the chord and oriented such 
that the jet exited in the spanwise direction. The jet blowing intensity was 
controlled by the jet supply plenum pressure and varied from 
C„ =0.0000 to 0.2236. 
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5.4.2 General Flow Features 

Figure 5.24 shows the close-up views of the velocity vectors at 10% chord 
downstream of the wing trailing edge for cases with spanwise wing-tip blowing 
off and on, and a = 8°. Figure 5.24 (a) shows clearly the formation of the 
trailing vortex at the wing-tip region. With moderately low blowing intensity, 
Figure 5.24 (b) shows that the trailing vortex is displaced not just outwards but 
also upwards. This is also observed in the experiment of Lee et al [1989]. 

Figure 5.25 shows a sequence of velocity vectors at the 90%-chord for 
ft = 2 with increasing blowing intensity. It is observed, especially in Figures 
5.25 (c) and (d), that in addition to the primary vortex, a secondary vortex 
rotating in the opposite direction of the primary vortex is formed at higher 
blowing intensities. The jet blowing becomes d ominan t and a typically free jet 
velocity profile is seen. This is also observed in the experiment of Lee et al 
[1989]. 


5.4.3 Forces and Moment 

Figure 5.26 shows the variations of the coefficients of lift, drag, pitching 
moment (about the wing quarter-chord) and LID with the blowing momentum 
coefficient for ft = 2 , 4 , 6 and 8°. It is observed in Figure 5.26 (a) that as the 
spanwise wing-tip blowing intensity increases, lift on the wing increases. The 
nonlinearity of the lift coefficient with the blowing momentum coefficient was 
also observed in the experiments (see Tavella et al [1985, 1988] and Lee et al 
[1989]). It is noticed that the computations predict the experimental results very 
well for cases of cr = 2 and 4 ; whereas, over-predictions by the computations 
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are observed for cases of a = 6° and 8\ Figure 5.27 shows the comparison of 
the computational results with the experimental data for the lift coefficient for 
various angles-of-attack without spanwise wing-tip blowing. Notice that no wall 
corrections to the experimental data were made. For low angle-of-attack cases, 
wall blockage is small and the computational and experimental results compare 
favorably. However, effects of the wall blockage become prominent at higher 
angle-of-attack and resulting in greater discrepancies between the computational 
and experimental results. 

Figure 5.26 (b) shows the computational results for the variation of the 
drag coefficient with blowing momentum coefficient. It is noticed that at 
moderately low blowing intensities, the drag increases. However, by increasing 
the blowing intensity further, the drag decreases. The drag coefficient may be 
lower than that without spanwise wing-tip blowing for high enough blowing 
intensity and lower angle-of-attack cases. It is also observed that for a particular 
blowing intensity, the drag increase is higher for higher angle-of-attack cases. 
This is clearly illustrated in Figure 5.26 (d) in which the variation of the LID 
ratio with the blowing intensity is presented. It is noticed that returns on LID via 
spanwise wing-tip blowing are more attractive for lower angle-of-attack cases 
than for higher angle-of-attack cases. 

Figure 5.26 (c) shows the variation of the pitching moment (about the 
wing quarter-chord) with blowing intensity. It is observed that as the blowing 
intensity increases, the pitching-up moment on the wing decreases marginally. 
Since spanwise wing-tip blowing is applied on almost the entire chord, the impact 
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of blowing on the pitching moment on the wing tends to even out and results in 
marginal changes in the pitching moment. 

5.4.4 Vortex Strength 

Figure 5.28 shows the variation of the streamwise distribution of the normalized 
total circulation with blowing momentum coefficient for the case of a = 4'. For 
streamwise stations immediately downstream of the wing trailing edge, the 
circulation integration loop basically encloses the tip vortex (or p rimar y vortex) 
and excludes the secondary vortex. Since spanwise wing-tip blowing has the 
effect of increasing the circulation strength of the primary and secondary 
vortices, a rather high normalized total circulation strength immediately 
downstream of the trailing edge is observed in Figure 5.28. However, as the 
vortices moves downstream, the circulation integration loop starts to enclose part 
of the secondary vortex which rotates in the opposite direction to that of the 
primary vortex. This results in a reduced circulation strength. For stations far 
downstream, the circulation integration loops eventually enclose both the 
primary and second vortices and an asymptotic value is reached for the vortex 
system. It is observed that the asymptotic circulation strength increases as the 
blowing intensity increases. This is intuitive since the overall lift on the wing 
increases as a result of spanwise wing-tip blowing, the total circulation strength 
for the vortex system far downstream must increase. 

Figure 5.29 shows the streamwise development of the pri mar y and 
secondary vortices. The formation of the primary and secondary vortices at 
high blowing intensity is clearly seen in Figure 5.29 (a). At the 90%-chord 
station, these vortices are distinctively apart. 
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5.5 Solution Accuracy 

Due to the constraints of the available computer speed and computer memory, 
computations were performed using 83x69x61 (i.e., 349,347) grid points. These 
constraints prevent a grid refinement study, which uses more grid points, to be 
carried out. However, the computational results obtained by using 51x51x41 
(i.e., 106,641) grid points are compared with the present results and the 
experimental data to study the effect of grid refinement on solution accuracy. 

Figure 5.30 shows the comparison in terms of C p distribution. It is 

observed that the leading edge suction peak can be better captured by the fine 
grid case as compared with the coarse grid case. The vortex induced low 
pressure peak near the trailing edge at the wing-tip region (see Figure 5.30 (b)) 
is also better captured by the fine grid case. Overall, the computational results 
obtained by the fine grid case have a better agreement with the experimental data 
as compared with those obtained by the coarse grid case. 

By comparing results obtained from the coarse and fine grid cases, a 
simple error analysis yields an error e(h) = [0(h)-0(ah)]/( 1-a 2 ), where (p 

denotes the solution and a denotes the grid refinement factor. Based on the lift 
coefficients obtained from the coarse and fine grid cases, the fine grid case gives 
an error of 2.5%. 
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(a) Side view. 
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(b) Plan view. 


Figure 5.2: Side and plan views of the trailing vortex 
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Figure 5.3: Velocity vectors at the quarter-chord showing flow from the lower 
surface to the upper surface of the wing near the wing-tip. 
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Figure 5.4: Distribution of the pressure coefficients, C p , on the wing near the 
wing-tip region for the case of a - 4°, = 0.17 and Re = 2.0 x 10 6 . 























Figure 5.6: Distribution of the pressure coefficients, C p , on the wing near the 
wing-tip region for the case of a = 12°, M_ = 0.17 and Re = 2.0 x 10 6 . 
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Figure 5.7: Formation of the tip vortex at x/c — 0.9 depicted by the stagnation 
pressure coefficient contours. 
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Figure 5.11. Variations of the coefficients of lift, drag and pitching moment 
(about the quarter-chord) for the case of Af_ = 0.17 and Re = 2.0 x 10 6 . 
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Figure 5.12: Streamwise distribution of the normalized total circulation for the 
case of a = 12°, M_ = 0.17 and Re = 2.0 x 10 6 . 
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Figure 5.13: Comparison of the calculated normalized downwash distribution 
downstream of the trailing edge with the theoretical analysis of Spreiter-Sacks 
for the case of a = 12°, M_ = 0.17 and Re = 2.0 x 10 6 . 
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Figure 5.14: The jet slot with a jet length of 11.2% chord and a jet width of 
1.5% chord is located at the 90%-chord location for the axial wing-tip blowing. 
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Figure 5.17: Variation of the streamwise distribution of the normalized total 
circulation with blowing momentum coefficient for the case of a = 8° 
M_ =0.17 and /te = 2.0xl0 6 . 
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(a) Without blowing. 



(b) With blowing. 


Figure 5.18: Particle traces showing the rolling-up of the trailing vortex with 
axial wing-tip blowing off and on. 
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Figure 5.20: Down wash distribution at various axial stations downstream of the 
tr ailin g edge for the case of cc — 8°, M_ = 0.17, Re = 2.0 x 10 6 and with axial 
wing-tip blowing on ( C M = 0.0047). 
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Figure 5.21. Variations of the coefficients of lift, drag and pitching moment 
(about the quarter-chord) with blowing momentum coefficient for the case of 
a = 8°, =0.17 and /te = 2.0xl0 6 . 



CHAPTER 5. RESULTS AND DISCUSSION 


96 



(a) Blowing off. 



(b) Blowing on ( = 0.0047). 


Figure 5.22: Pressure contours at x/c = 0.95 near the wing-tip region for cases 
with axial wing-tip blowing off and on. 
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Figure 5.23: Wind-tunnel model used in the experimental studies with spanwise 
wing-tip blowing (taken from Tavella et al [1988]). 
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(a) Blowing off. 





(b) Blowing on ( = 0.0169). 

Figure 5.24: Close-up views of the velocity vectors at x/c = l.l for the case of 
a = 8° with spanwise wing-tip blowing turned off and on. 
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(a) C M =0.0000. 



(b) C, =0.0168. 


Figure 5.25: Velocity vectors at xfc = 0.9 for the case of a = 2* with increasing 
span wise wing-tip blowing intensity. 
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(c) C M =0.0342. 
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(d) C; =0.2076. 
Figure 5.25: Concluded. 
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Figure 5.26: Variations of the coefficients of lift, drag, pitching moment (about 
the wing quarter-chord) and LID with blowing momentum coefficient for 
a = 2\ 4\ 6° and 8*. 
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(c) Pitching moment coefficient (about the wing quarter-chord). 



Figure 5.26: Concluded. 





Lift Coefficient 





2 r/cv_ 



Figure 5.28: Variation of the streamwise distribution of the normalized total 
circulation with blowing momentum coefficient for the case of a- 4\ 
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(a) jc/c = 0.95. 



(b) x/c = 1.00. 


Figure 5.29: Stagnation pressure contours showing the streamwise development 
of the primary and secondary vortices for the case a- 2* and C„ = 0.2067. 
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(c) xlc = 2.00. 



(d) xlc = 3.00. 


Figure 5.29: Continued. 
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(a) At 2ylb = 0.676. 



(b) At 2y/b = 0.971. 


Figure 5.30: Comparison of the computational results using coarse and fine 
grids with the experimental data for the case of a = 12*. 






Chapter 6 


Conclusions and 
Recommendations 


6.1 Conclusions 

N ume rical simulations of the complex flowfield of the trailing vortex of a wing 
with axial and spanwise wing-tip blowing have been successfully performed 
using Computational Fluid Dynamics techniques. The unsteady, three- 
dimensional, thin-layer Navier-Stokes equations are solved using a time-accurate, 
im plicit, finite difference numerical algorithm. The actuator plane concept and 
the solution blanking feature of the chimera scheme have been successfully 
incorporated to simulate wing-tip blowing. Computational results are generally 
in good agreement with experimental data. The main conclusions of the present 
study are summarized as follows: 

(1) The axial wing-tip blowing has the effect of modifying the 
flowfield behavior immediately downstream of the wing. 
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However, its effect diminishes farther downstream. Overall, the 
axial wing-tip blowing strengthens the circulation strength of the 
trailing vortex marginally due to the marginal increase in the wing 
lift. 

(2) At moderately low blowing intensity, the spanwise wing-tip 
blowing has the effect of displacing the tip vortex outboard and 
upwards and results in the augmentation of lift. 

(3) At higher blowing intensity, a secondary vortex rotating in the 
opposite direction to the primary vortex is formed by the spanwise 
wing-tip blowing. 

(4) The circulation strength of the vortex system increases as a result 
of lift augmentation on the wing due to the spanwise wing-tip 
blowing. 

(5) Spanwise wing-tip blowing can be an effective means of control to 
effect a change of vortex position on one wing-tip and thereby 
produce a rolling moment through asymmetric wing lift. 

6.2 Recommendations 

Some recommendations for future work are as follows: 

(1) Improve the solution accuracy by increasing the grid point density 
in the regions of the leading edge, trailing edge, tip of the wing and 
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the trajectory of the trailing vortex. Study the effect of grid 
refinement on solution accuracy. 

(2) Extend the outflow boundary to several wing spans to capture the 
far-field behavior of the trailing vortex. 

(3) Perform full Navier- Stokes simulations to investigate the 
limitations of the thin-layer approximation for the trailing vortex 
flow. 

(4) Model the complex flowfield of the trailing vortex with wing-tip 
blowing with better turbulence models. 

(5) Simulate the complex flowfield of the trailing vortex with wing-tip 
blowing for a rotorcraft by incorporating the rotating motion of 
the rotor blades. Fine grid points must be clustered to capture the 
trailing vortices from the rotor blades as they spiral downwards 
due to the rotor downwash. 
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